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Studies of large systems of interacting oscillatory elements is a popular and 
extensively developing branch of nonlinear science. The number of publica¬ 
tions on the snbject grows rapidly, with many crucial contributions pnblished in 
the Chaos journal. In addition to pnrely academic interest, this research finds 
promising applications in varions fields. Representative examples are nnder- 
standing of pedestrian synchrony on footbridges and of other social phenomena, 
development of efficient high-frequency power sonrces, modeling and control of 
nenronal rhythms, etc. In this paper we present our view of the recent devel¬ 
opment of this highly interesting field. 


I. INTRODUCTION 

In the second half of the 17th century Engelbert Kaempfer, Dutch physician, made a jour¬ 
ney to South-East Asia and later published a book, describing his trip>i. In particular, in 
his memoirs he gives an account of a spectacular show, which happens if a swarm of fireflies 
occupies a tree: the insects “hide their Lights all at once, and a moment after make it appear 
again with the utmost regularity and exactness”. This phenomenon of self-synchronization 
in a large population of interacting oscillatory objects not only remains an appealing enter¬ 
tainment — be it an excursion on a night river in Thailand to observe fireflies or cycling 
in a large group of people, equipped with electronic “bikefiies”, as a part of a festival in 
Chicago — but it stays in the focus of scientific interest within many decades. The studies 
of various aspects of the collective dynamics in large oscillatory networks attract attention of 
physicists and applied mathematicians and find applications ranging from electrochemistry 
to quantum electronics, and from bridge engineering to social sciences. 

A particularly popular and expanding field of applications is neuroscience. For the first 
time the link between synchronous hashing of the fireflies and origin of the brain rhythms 
was established, on a quantitative level, by the famous mathematician Norbert Wiener in 
his monograph on Cybernetics^, in the chapter entitled “Brain waves and self-organizing 
systems”, see also Ref. j^. Putting forward the hypothesis that the brain waves emerge due 
to “phenomenon of the pulling together of frequencies”, he questioned, whether the same 
nonlinear mechanism takes place in case of fireflies, crickets, and other species exhibiting 
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collective oscillatory behavior, and suggested that experiments on hreflies and on electronic 
systems can shed light on the brain wave dynamics. 

Since that many experiments have been reported, including those with electrochemical^*^ 
and electronic oscillators^, metronomes^, Josephson junctions^, laser arrays^, yeast cells^^, 
and gene-manipulated clocks in bacteriap^. They are complimented by observations of many 
social phenomena like pedestrian synchrony on footbridges^, synchronous hand clapping in 
opera houses^, egg-laying in bird colonies^, and menstrual synchronyi^ (the latter effect 
remains controversial^). This research was also accompanied by an essential progress in the 
theoretical description, that we outline below, along with the open questions. 


II. SOLVABLE MODELS 


A few years after publication of the Wiener’s book, Arthur Winfree presented a hrst math¬ 
ematical description of collective synchrony in a large population of biological oscillators^. 
Reducing the dynamics of each oscillator to that of only one variable, the phase (we will 
discuss conditions for this reduction below), he proposed the mean-held model 

N 

( 1 ) 


= Wfc + ^r((pfc) ^ J((pj), 
i=i 


where iV ^ 1 is the number of units, ujk are their natural frequencies, and e quantihes 
the strength of the interaction. The function r(<yCfc) describes the phase sensitivity of the 
oscillator to an inhnitesimal perturbation and is typically called the phase response curve, 
PRC. Notice that the PRC can be experimentally obtained by repeating stimulation of an 
isolated system. Finally, the forcing function describes the effect of the j-th unit on 

the unit k. In this setup the coupling is assumed to be global, i.e. of the all-to-all type, 
and functions P, I are assumed to be identical for all interactions. Thus, the inhomogeneity 
of the system is due to a distribution of frequencies only. Although the model is quite 
complicated for the analysis, Winfree has shown that it exhibits a transition to a macroscopic 
synchronized state, characterized by non-zero mean held N~^ 'Yhj He discovered that 

the collective synchrony is a threshold phenomenon; the transition occurs if the coupling 
strength e is large enough or the inhomogeneity, i.e. the width of the distribution of Uk, is 
sufficiently small. For recent studies on the Winfree model see Refs. it has been 


shown that it can be treated analytically at least if functions F, I contain only one Fourier 
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harmonics. 

The next pioneering step has been done by Yoshiki Kuramoto in his seminal publication 
40 years agc3. The model he suggested can be considered as the weak-coupling limit of 
Eq. ([T]). Indeed, if e <C cu, then each Eq. ([T]) can be averaged over the oscillation period; 
then each term yields a function of the phase difference, gi’-Pj — ’-Pk)-, so that the 

averaged model reads 

N 

- (Pk) ■ ( 2 ) 

i=i 

The general case of an arbitrary 27r-periodic function g is discussed later, while the simplest 
case g{-) = sin(-) corresponds to the famous Kuramoto model (notice that choice of the sine 
function is not a result of some approximation but just the simplest solvable case): 

N 

= cufc + — ^ sin(y?j - Lpk) =uJk + eRsm{Q - ipk) • (3) 

t=i 

Here RR® = N~^ is the complex mean held, R and 0 are its amplitude and phase, 

respectively. Kuramoto solved the problem in the thermodynamical limit N ^ oo, using 
a self-consistent approach: assuming a harmonic mean held with unknown amplitude and 
frequency, he obtained closed integral equations for these two quantities. The celebrated 
result is the existence of the critical coupling, Ec, proportional to the width of the frequency 
distribution. For sub-threshold coupling the mean held is exactly zero, while for e > 
a nontrivial solution with non-zero mean held appears; the amplitude of the held grows 
as \/£ — £c and its frequency equals the central frequency of the distribution of Uk (which 
is assumed to be symmetric and unimodal). Thus, appearance of the collective mode can 
be treated as a second-order nonequilibrium phase transition. As has been shown later, 
the Kuramoto model with the uniform frequency distribution exhibits a jump of the order 
parameter—. 

III. COLLECTIVE DYNAMICS OF THE KURAMOTO MODEL 

The results of Kuramoto gave an enormous impact on the development of the held, with 
still an increasing number of publications on the subject. The model ([3]) and its extension 
due to Sakaguchi and Kuramoto^S, who accounted for a possible phase shift in coupling, 

= cufc-h £Rsin(0 - -F a) , (4) 
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became a paradigmatic model for analysis of large oscillator ensembles. 


A. Watanabe-Strogatz theory 


We now briefly discuss an interesting and important mathematical property of the 
Kuramoto-Sakaguchi model (jlj), namely its partial integrability. We start this discussion 
by consideration of the simplest case of identical units, ook = oj- As has been shown in the 
seminal publications by Watanabe and Strogatz^^ (WS), the dynamics of this system can 
be described by three global variables p, $, and N — 3 constants of motion ipk. Here the 
variable 0 < p < 1 is roughly similar to the mean held amplitude R] <h and ^ are angular 
variables; often it is convenient to combine two variables as z = pe‘‘*’, see also Refs. 241. 
The main idea of the powerful WS theory is as follows. Consider N > 3 identical oscillators 
governed by 

(Pfc = uj{t) + Im [H , 


where H(t) is an arbitrary common forcing. Obviously, Eq. 
The latter can be re-written as 


(5) 

is a particular case of Eq. (jSD. 


d 1 

- (e'*’*) = + -Hit) - —. 

Next, we transform N variables (fk to complex z, \z\ < 1, and N real ^k, according to 

1 -|- z*e*5'= 

This transformation, found by WS and written in this form in Refs. 


( 6 ) 


(7) 


24 


25| , is known as the 


Mobius transformation. Since the system is over-determined, one requires 

N 

_ /okfc\ _ 


7V-i^e^^'= = (e^^'=) = 0 . 


( 8 ) 


k=l 


Substituting Eq. ((^ into Eq. 


i + 


we obtain, after straightforward manipulations: 

H 


zz* — zz* + i'Cfc(l — l^r) 


e'^'' = icjz + 


H* 


- —z^ + [ia;(l + l^n + {z*H - zW 


2 

}ik 


(9) 


+ 


H 


R -f icaz" -F —2;*^ - 


H* 
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Averaging these equations over k, using Eq. 

H H* 


Z — \(jJZ — 


-z^ = 


2 2 

and this equation is obviously satished, if 


and = 0, we obtain 

H* H 




H H* 


*2 


(e2'«'=) 


i = icj^; + „ - 

2 2 

Substitution of Eq. ffTOj) and its complex conjugate into Eq. 

z*H - zH 


( 10 ) 


i&(i-|zr) = 


lU + 


yields: 


Excluding the fully synchronous case \z\ = 1 and introducing = « + t/'fc, where all 
'ifjk = const, we obtain: 

a = oj + \m.{z*H) . (11) 

Expressions fllOlllip represent the Watanabe-Strogatz equations, which completely describe 
the evolution of an ensemble of identical oscillators. 

Extension of the WS theory for the case of non-identical units depends on the structure of 
the ensemble. Consider hrst the case of a hierarchical population, which consists of M groups 
of identical units, so that units in each group are subject to the same held^. In this case the 
dynamics of the ensemble obeys the system of M coupled WS equations fllOlllip . Another 
practically important case is a large population, N ^ oo, which can be characterized by a 
distribution of frequencies g{u)). This system is described by WS variables z{u),t), a{uj,t) 
and the equations for dtz, dta are similar to (I10|llip . see Ref. [24 1. 


B. From WS to Ott-Antonsen theory 


There exist a particular case, when the WS equations can be essentially simplihed. As 


has been shown in Refs. 


24j |. for large N and for the uniform distribution of constants of 


motion ijjk, the WS variable coincides with the local Kuramoto mean held 


i>2tt 


Z{u,t) = / w{(p,t\u)e^'^d(p , 


( 12 ) 


(where w{(p, t\u) is the distribution density of oscillators’ phases at frequency u] f dipw{(p, t\u) 
1) to be distinguished from the global mean held 


^{t)= / g{uj)Z{cj,t)dip . 


(13) 
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( 14 ) 


As a result, one obtains closed equations for Z\ 

-- + — - 

+ . (15) 

In the most common case of the mean held coupling H = ee^^Y, the forcing H is independent 
on a and Eq. (IT^ becomes irrelevant. Hence, we are left with Eq. (fT4|) which also appears 
in the recent theory by Ott and Antonsen (OA)^>21^ briehy introduced below. 

Consider the thermodynamical limit of the model 




(16) 


and the corresponding continuity equation for the distribution of phases 

dw d 

_ + _(„^) = o. 


(17) 


Next, let us introduce Fourier components of the density (the generalized local Daido order 
parameters) 

p27V 

Zm{uj,t)= / w{uj,ip,t)e""^'^dLp , (18) 


cf. Eq. flT^ . Computing 


dZ{u), t) 
dt 


f27r 


dw : 


f27r 


dt 


e™^d(/3 = - 


d{WLp) ; 


dt 


e'"*^d(/3 , 


integrating by parts and using Eq. (fT6D . one obtains an inhnite-dimensional system of ODEs: 

(19) 


dZ{uj,t) _ . m 

yWjOjZYn ~\~ 2 \HZyyi—\ H ZyjiJ^ij 


A particular case Zm = = Z"^, called the OA manifold, reduces all the equations (IT^ to a 

single Eq. flT^ . Thus, the OA manifold corresponds to the special solution of the WS theory, 
with the uniform distribution of constants of motion ij). Furthermore, OA argued that for 
continuous frequency distribution g{uj) the OA manifold is the only attractor— (although 
relaxation to it may be rather slow^^). A particular case of the Lorentzian distribution 
g{u) = [7r(a;^ + 1)]“^ admits a further essential simplihcation. Under assumption that Z{u) 
is analytic in the upper half-plane, one computes the integral in Eq. flT^ by virtue of residues 
and obtains Y = Z(i). Substituting this along with uj = i into Eq. (1T4)) one obtains the OA 
equation 


Y 



-1 F- 


ee 


-i/3 


-Y^Y* 


( 20 ) 
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for the evolution of the mean held in the Kuramoto-Sakaguchi model. Looking for a station¬ 
ary synchronous solution, we set Y = and obtain the amplitude Rq = y/l — 2/e cos {3 

and frequency v = [e cos /5 — 1) tan 13 of the mean held. 

For an illustration of the WS and OA theories we consider the model of two interacting 
populations by Abrams et ah—. Both populations consist of the same number of identi¬ 
cal oscillators, but the coupling strength within the group dihers from that between the 
groups. The dynamics can be fully described by two coupled systems of WS equations, 
and is therefore six-dimensional2i. In the particular case of uniformly distributed constants 
'ijj, i.e. on the OA manifold, the dimension reduces to four. The latter case, studied in 
Ref. 29j, reveals, for certain parameters, an interesting solution, when one population syn¬ 
chronizes while the other stays between complete synchrony and full asynchrony, i.e. = 1 
and 0 < i ?2 < 1- This symmetry-breaking state is called chimera. (Notice that originally 
chimera states have been introduced for nonlocally coupled oscillator chains22.) Interesting, 
that i ?2 can be time-periodic. Analysis of the full six-dimensional systems exhibits addi¬ 
tional solutions with the quasiperiodic chimera states^i. Theoretical predictions have been 
conhrmed in recent experiments with two groups of metronomes, placed on platforms which 
were coupled via springs^. 


C. Generalizations of the Kuramoto model 


There are many studies of different generalizations of the Kuramoto model. Here we 
briefly mention those where the coupling function is purely harmonic like in Eq. (j4]), but the 
overall setup is more complex: 

a. Multi-modal frequency distribution and several interacting populations. We have 
seen that an ensemble with a Lorentzian distribution of frequencies is described by the 
OA Eq. fl20|) . A multi-modal distribution of frequencies can be modeled as a superposi¬ 
tion of Lorentzians and, hence, described by a system of coupled OA Eqs. (l20ll . see e.g. 
Ref. {mI- Moreover, this approach can be generalized to a set of populations with frequen¬ 
cies, distributed around completely different central values, whereas the latter can be either 
in resonance^ or not^. 


b. Nontrivial transitions for unimodal distributions. For a long time it has been as¬ 
sumed that in the Kuramoto-Sakaguchi model (jl]) with different unimodal distributions of 





frequencies the dynamics of the mean held is qualitatively the same as for the Lorentzian 
distribution (Eq. fl20|) h Rather surprisingly, Omelchenko and Wolfrum^i have demonstrated 
rather complex transition scenaria, including hrst-order transitions and bistability, for some 
unimodal distributions. 


c. Complex coupling schemes. The Kuramoto-Sakaguchi model describes a “direct” 
coupling scheme: the mean held, calculated algebraically from the states of all oscillators, 
enters the equations for the phases. The coupling scheme can be, however, more complex: 
the mean held may act on some macroscopic variables that obey a set of generally nonlinear 
diherential equations, and the acting force is a function of these variables. For example, 
in a description of pedestrian synchrony on a footbridge^, one describes each pedestrian 
by an individual phase variable, but one needs also equations for the swinging mode of the 
bridge. The latter is driven by the held created by all pedestrians and, in its turn, ahects 
their gaits. Similarly, electronic)^ or electrochemical^i^ oscillators can be coupled through the 
common macroscopic current or voltage, which obeys macroscopic equations describing the 
coupling circuit. In this way one also describes synchronization of Josephson junctions^ or 
spin-torque^ oscillators. 


d. Nonhomogeneous populations. There is another generalization of the standard 
Kuramoto-Sakaguchi coupling. The latter assumes that all the oscillators make the same 
contribution to the mean held and that the mean held acts on all oscillators in an equal 
way. Refs. 


38 and 


39l | generalized this to the situation, where the global held still can be 


introduced, but oscillators contribute to it diherently, i.e. with diherent amplitudes and 
phase shifts, and the held also acts diherently on diherent oscillators. For a physical imple¬ 
mentation one can can consider a receiver which collects signals emitted from the oscillators 
(where the attenuation and the phase shift are due to signal propagation properties), and 
the governing signal is then transmitted to the oscillators, being also subject to attenuation 
and phase shift depending on the positions of the units^. One variant of an inhomogeneity 
of the population is when it consists of two groups with diherent reaction to the force^: 
some are “conformists” (they follow the force) and some are “contrarians” (they tend to be 
in anti-phase with the forcing). Another possibility, quite general for neural ensembles, is 
when some oscillators are inhibitors (trying to desynchronize others, coupling is repulsive) 
and others exert an excitatory action, attempting to synchronize (attracting the phase )^‘^. 
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e. Effects of noise. Independent noisy forces acting on oscillators of a population coun¬ 
teract synchronization. In this sense noise is a source of disorder, similar to the distribution 
of natural frequencies. Due to noise, synchronization is a threshold phenomenon already 
for identical oscillators and the transition occurs at a critical coupling which is proportional 
to the noise intensity. In the thermodynamic limit such a system is described analytically 
by a nonlinear Fokker-Planck equation which differs from the continuity equation (fT7|) by a 
term ~ where a is the amplitude of the noise. Completely opposite is the action of 

a common noise: it tends to synchronize the population of oscillators and for identical units 
this effect can be described within the WS framework^. 

/. Finite-size fluctuations. Kuramoto and OA theories have been developed in the ther¬ 
modynamic limit of inhnitely large populations; WS theory applies to any population size, 
but for identical populations only. In hnite populations with different natural frequencies 
of units one expects to observe fluctuations, both prior and beyond the synchronization 


transition, which is dehned ambiguously in this case. In Ref. 


44| the Kuramoto model with 


a uniform distribution of frequencies and a relatively small number of oscillators have been 
shown to be chaotic prior to synchronization transition, the maximal Lyapunov exponent 
however decreases with the system size as A ~ N~^. Above synchronization transition only 
regular dynamics have been observed. However, for A^ S> 1 and close to the synchroniza¬ 
tion transitions the regime is complex: if it is not chaotic, then it is quasiperiodic with a 
large number of incommensurate frequencies; here statistical approaches based on hnite-size 
scaling have been applied to hnd critical indices of A^-dependence^*^. 

g. Kuramoto model on networks. Kuramoto model has been initially formulated for the 
ensemble of globally coupled oscillators. Recently, it has been extensively studied for other 
coupling conhgurations, i.e. for networks of different complexity, including small-world net¬ 
works^, multidimensional hypercubic lattices^, networks with a modular structure^, and 
an ensemble with an extra leading element (hub)^. Dynamical properties of the transition 
to synchronization depend on the network topology^. One of the popular applications here 
is study of synchronization of power grids^. 

h. External forcing and collective phase resetting. As one can see from Eq. (120]) . the 
macroscopic order parameter obeys an equation for a self-sustained oscillator. Thus, the 
collective mode has the same properties as such an oscillator. In particular, if the ensemble 
is forced by a periodic force, the latter can entrain the frequency of the mean held oscilla- 
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In the case of a pulse force, the latter can shift the phase of the collective mode; 
in this way one defines the phase response curve for the collective mode^. 

i. Mathematical results. The Kuramoto model has been extensively studied on the 
physical and computational level, but rigorous mathematical results are sparse. They mainly 
refer to stability analysis of the desynchronized state2^>^; for a description of the synchro¬ 
nization transition as a bifurcation problem see Refs. j^ |. 


IV. GENERAL COUPLING FUNCTIONS 


Now we come back to Eq. ([2]). For a general coupling function g, this equation represents 
the Daido model^'^. Expanding g into Fourier series, g{r]) = ^ introducing 

generalized order parameters 

N 

( 21 ) 

i=i 

one can re-write the model as 






( 22 ) 


One can see that generally all order parameters should be determined self-consistently, so 
that a complete analysis at general coupling is still missing. We mention here several inter¬ 
esting regimes appearing due to high harmonics in coupling function. 

j. Clustering. Even for identical oscillators the WS theory does not apply, and a pop¬ 
ulation can build several clusters^, each of them consisting of fully synchronized units. 

k. Heteroclinic cycles. In Ref. j^] nontrivial regimes of clustering and switching be¬ 
tween different cluster states have been found in a population of identical oscillators with a 
function g containing the first and the second harmonics. This complex dynamics is due to 
a heteroclinic cycle in the phase space, well understood for small networks^. 

l. Multi-branch entrainment. Already for two harmonic components in the coupling 
function g, the r.h.s. of Eq. (l22l) as a function of can have two stable branches. This 
is a new situation compared to the standard Kuramoto model: now for a given mean field 
entrainment at two different microscopic phases is possible. This leads to an enormous 
multiplicity of microscopic states and to a complex structure of macroscopic regimes^'^'^. 
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V. NON-PHASE MODELS 


We started our discussion of collective ensemble dynamics from phase models. This 
approach relies on the well-known idea that motion along the limit cycle of an autonomous 
system can be parameterized by a single variable, the phase. Moreover, if the interaction 
of the oscillator with the environment is weak so that one can neglect the deviation of 
the trajectory from the cycle of the autonomous system, then the low-dimensional phase 
description remains valid. Generally, for strong coupling, one has to analyze full dynamical 
models which is a complicated problem that can hardly be treated analytically. 

Numerical studies reveal that transition to collective synchrony is a quite general property, 
observed for various periodic, noisy, and even chaotic oscillators^!^, including, e.g., spiking 
and bursting model neurons. The picture is, however, non-universal. The most transparent 
and studied model is ensemble of coupled Stuart-Landau oscillators^ (this oscillator is the 
simplest prototype of a limit-cycle oscillator, with a perfect separation of amplitude and 
phase variables, see Eq. (12T|) below), and essentially new effects are the oscillation quenching, 
when too strong coupling effectively introduces additional damping to ensemble elements, 
and the collective chaos in a system of initially periodic units. 

On the other hand, non-identical chaotic phase-coherent Rossler oscillators adjust their 
frequencies and phases (this effect is known as phase synchronization of chaos^) and pro¬ 
duce a nearly periodic mean held. The oscillators themselves remain chaotic, but irregular 
huctuations of the amplitudes turn to be averaged out in the mean held; the small huctu- 
ations of the latter are presumably due to the hnite-size ehect^. Experimental studies on 
chaotic electrochemical oscillators^ conhrm theoretical predictions. 

Another important class are rotators: these systems are described by an angle-like vari¬ 
able, which is very similar to the phase, and do not have amplitudes. If an equation governing 
the rotator’s angle is one-dimensional, the dynamics can be reduced to a Kuramoto-type 
phase model, what has been extensively discussed in context of Josephson junctions^^. Quite 
often the inertia of rotators cannot be neglected and, hence, they are described by a second- 
order equation for the angle variable. In this case the dynamics can strongly deviate from 
that of the Kuramoto model, e.g. transition to synchrony may exhibit hysteresis, similarly 
to a first-order phase transition^. A particular subclass is constituted by globally cou¬ 
pled rotators without damping: this conservative system yields the so-called Hamiltonian 
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68l | for a review of recent results and relation to the 


Mean Field (HMF) model^; see Refs. 

Kuramoto model. 

VI. COMPLEX COLLECTIVE DYNAMICS AROUND SYNCHRONY 

We have discussed in detail the main effect observed in globally coupled systems, namely 
emergence of the collective mode, which is well-understood in the simplest case, when many 
units synchronize and therefore there outputs sum up coherently, contributing to this mode. 
We have also mentioned, that the mode itself can exhibit chaotic dynamics. Now we discuss 
other, less explored, situations, when the collective dynamics of an ensemble is complex. 

A. Partial synchrony 

In case of the Kuramoto model of identical oscillators, only two situations are possible: 
full synchrony for attractive coupling (order parameter equal to one) and fully asynchronous 
state (zero order parameter) for repulsive coupling. However, this situation is not general 
and we can face a case, when both fully synchronous and fully asynchronous states are 
unstable, so that the system settles somewhere in between, at a state which is often called 
partial synchrony. The simplest and well-known example of partial synchrony is clustering; 
below we discuss several other situations where oscillators are not organized in clusters. 

We examplify partial synchrony with iV S> 1 mean held coupled oscillators^: 

Xk = yk-xl + 3xl- Zk + 5 + e{X - Xk) , 

yk = 1-5x1-yk , (23) 

Zk = 0.006[4(xfc -I- 1.56) - Zk] , 

where k = and X = Here individual units represent a popular 

Hindmarsh-Rose neuronal model^. The quite standard parameter values used here corre¬ 
spond to a limit-cycle solution for the uncoupled neurons, or, in neuroscience language, to a 
state of periodic spiking. Fully synchronous state, Xi = X 2 = ■ ■ ■ = x^, 2/i = 2/2 = • • • = yN, 
zi = Z 2 = ■■■ = zn is obviously a solution of the system. However, stability of this state 
depends on the coupling strength £, as can be checked numerically by virtue of computa¬ 
tion of the transversal stability via evaporation multipliers^!^. This analysis, conhrmed by 
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direct simulation, demonstrates, that there exist a critical coupling value at which the 
synchronous solution loses its stability via a Hopf-like bifurcation, i.e. two complex mul¬ 
tipliers cross the unit circle, giving birth to a new frequency. (Notice that the stability of 
the synchronous solution is re-gained for very large e.) Beyond the synchrony breaking, the 
states of oscillators in the phase space form a thin stripe, stretched along the limit cycle; the 
points within this stripe slowly interchange their position, with a characteristic time of tens 
of periods. Roughly speaking, the average frequency of all oscillators and of the mean held is 
the same, but the phase shift of the oscillators with respect to the held is slowly modulated. 
As a result, the dynamics looks quite complicated and is possibly weakly chaotio^. 

For another example we consider a popular model of a series array of resistively shunted 
Josephson junctions. The junctions are coupled by virtue of a parallel RLC-loa.d^. As has 
been shown in cited refs., for a weak coupling and linear load, the system is equivalent to 
the Kuramoto model. Consider now a nonlinear coupling] namely, let the inductance in the 
RLC-circuit be nonlinear so that the magnetic hux depends on the current Q through the 
RLC-load as $ = LqQ + LiQ^. Numerical study^ reveals a synchrony-breaking transition; 
at e = Ec the synchronous state becomes unstable; at this critical coupling value real evap¬ 
oration multiplier p becomes larger than one (in contradistinction to example fl2^ where 
the multiplier is complex). For e > Sc the systems is in a state of partial synchrony, with 
the order parameter being a smooth decreasing function of e. Furthermore, the dynamics 
becomes quasiperiodic: the frequency of the mean held is larger than the frequency of in¬ 
dividual junctions, so that the latter are not locked to the held. The frequency diherence 
grows with e — Ec- 

A general theory of partially synchronous states which appear after the synchrony break¬ 
ing is still missing and requires further investigations. Below we present an analytically 
tractable and relatively transparent example. 

B. Self-organized quasiperiodic dynamics 

Now we analyze the system of iV ^ 1 Stuart-Landau oscillators: 

hfc = (1 -l- iuj)ak — (1 -|- ia)\ak\‘^ak (^i -l- i^2)R (24) 

-(hi +iV2)\A\‘^A , 
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where A = N ^ '^f=i .^ 1^25 hi ,2 are coupling parameters. This model differs from that 


of Refs. 6 J] due to the nonlinearity in coupling. In the weak coupling approximation the 
model with purely linear coupling^ (i.e. with hi = h 2 = 0) yields the Kuramoto-Sakaguchi 
Eq. dlj), while the nonlinear Eq. fl2T|) reduces to a particular case of the following phase 
model^: 


(p = u + ea{e, R) sin[0 - pk + (3{e, R)] 


(25) 


Equation can be considered as an extension of the Kuramoto-Sakaguchi model. Here R 
is the mean held amplitude and the bifurcation parameter e corresponds to a combination of 
the parameters ^ 1 , 2 , hi ,2 and a{e, R), R) are some functions. Notice that Eq. fl25|) appears 
also in a model of Stuart-Landau oscillators coupled via a common nonlinear medium^>^, 
similarly to the Josephson junction model. 

Let us consider the effect of these functions separately, starting with the case when 
(3 = const, |/3| < n/2, i.e. the coupling is attractive, cf. Refs. HI Q. Suppose that a is 
a decreasing function of e, e.g. a{e, R) = {1 — eR^)R (this function indeed appears for a 
certain combination of parameters in Eq. (I24l) ). For e < Ec = 1 this function is positive for 
fully synchronous case R = 1 and, hence, this state is stable. For e > EcWe have a{e, 1) < 0, 
i.e. the coupling becomes repulsive. As a result, the system stays at the border between 
attraction and repulsion, determined by the condition sR^ = 1, forming a self-organized 
bunch stated. In this state, for general initial conditions the oscillator phases spread around 
the unit circle so that R = this bunch is stationary in the coordinate frame, rotating 

with the frequency u. 

Now we analyze a more interesting case when a = R and {3{e,R) = (3o + |/9o| < 

7r/2, /9i > 0. It is easy to see that synchrony (i? = 1) is stable if /3o+/Si£^ < tt/2 and becomes 
unstable when e exceeds the critical value Ec = -s/ /2, — /9o)//^i- Again, the system settles 
at the border of stability, so that the condition /Sq + jdiE'^R^ = 7r/2 is fulhlled. However, in 
this case the dynamics exhibit a peculiar feature, also possessed by the Josephson junction 
model discussed above: the frequency of the mean held differs from the frequency of the 
individual units. Thus, the state can be characterized by two generally irrationally related 
frequencies, and therefore is denoted as self-organized quasiperiodicity (SOQ). Qualitatively, 
the emergence of quasiperiodic motion can be explained as follows. For e > e^ the system 
is partially synchronous, i.e. R < 1 and all oscillators have different phases (notice that for 
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general initial conditions without clusters the phases must be different according to Eq. 0). 
Hence, the instantaneous frequencies of oscillators differ and, therefore a stationary (in a 
rotating coordinate frame) distribution is not a solution. Quantitative analysis of SOQ 
dynamics with computation of the frequencies of the collective mode and of the oscillators 
can be found in Refs.— SOQ states in real-world systems have been demonstrated in 
experiments with electronic circuits with global nonlinear coupling^. 

To complete the discussion of this issue we mention that similar complex states with a 
nonzero collective mode can appear also without desynchronization transition. For some 
systems, the fully synchronous and completely asynchronous states are unstable already for 
inhnitely small coupling. A well-studied example is the van Vreeswijk model of globally 
coupled leaky integrate-and-fire neurons^. Another example is the emergence of dephasing 
and bursting in a system of Morris-Lecar neuronal models^; computation of the evaporation 
multipliers for this system shows that the synchronous state is unstable for the positive 
coupling range, where the effect is observed. 

C. Chimeralike states in globally coupled systems 

We have already mentioned a symmetry-breaking chimera state in a system of two inter¬ 
acting Kuramoto populations. Now we discuss emergence of a chimeralike state in a single 
homogeneous population. At first glance, such regimes are not allowed, because identical 
units subject to a common force should exhibit the same dynamics (i.e. to be either all 
synchronized or all desynchronized). On the other hand, there is a number of numerical 
observations of the mixed states, where only a fraction of the population merges to one or 
several clusters, while other elements remain scattered^, see also recent studies of chimera 
states in linearly and nonlinearly coupled Stuart-Landau oscillators^. The conditions for 
emergence of such mixed states are not yet clear. Nevertheless, we can outline one mecha¬ 
nism which results in splitting of the population into coherent and incoherent groups. 

Identical elements subject to the same force can behave differently if they are bistable 
(multistable), i.e. possess at least two nontrivial attractors. Another requirement is that 
the mean field forcing shall be synchronizing for the units on one attractor and repulsive for 
the elements on the other one. Taking into account that this partially coherent and par¬ 
tially incoherent state shall be maintained self-consistently, we conclude that this condition 
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is not trivial. To illustrate this mechanism, we hrst consider a rather artificial but trans¬ 
parent example^, where the oscillators are the non-isochronous modified Stuart-Landau 
systems. The modification refers to the nonlinearity: in addition to the 3rd order term 
lo/cpOfe we add also the terms ~ \ak\‘^ak, ~ \ak\^ak, cf. Eq. fl2T|) . As a results the systems 
possess two stable limit cycles with d ifferent frequencies, ^2 > The global coupling 
has its own dynamics, cf. Refs. 
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80| . so that the mean field forces a harmonic oscillator. 


U + JU + rfu = N~^ Re(aj), and its output ii acts on the Stuart-Landau systems. Sup¬ 
pose now that parameters are chosen in such a way that hli < rj < fl2- Since the phase shift 
introduced by the harmonic oscillator in the last equation is frequency-dependent, the cou¬ 
pling synchronizes the large amplitude limit-cycle oscillations, but prevents synchronization 
of the low-amplitude ones. 


Much less trivial is emergence of coherent-incoherent states in an ensemble of globally 
coupled oscillators with internal delayed feedback loop. Such oscillators naturally appear, 
e.g. in laser optics^ as well as in numerous biological applications^. In the simplest case 
the autonomous dynamics of an oscillator with a delayed feedback loop can be described by 
a phase model, ip = u + asin((/?.r — 9?), where ipr = p>{t — r), r is the delay, and a quantifies 
the feedback strength^i^. Assuming the global coupling in the ensemble of such oscillators 
to be of the Kuramoto-Sakaguchi type, one writes the model as^: 


= a; + a sm{ipr,k - Vk) + eRsm{Q - + f^) ■ (26) 

Stability analysis of the fully synchronous one-cluster state yields that it is unstable for 
(3 > 7r/2. However, numerical simulation shows that for f3 > 7r/2 the asynchronous state 
with zero mean field, i? = 0, is also unstable (see^ for other parameter values). Thus, 
the system “chooses” a state between full coherence and full incoherence, and in a range 
of parameters this state is reminiscent of a chimera: there is one big cluster (about 70% of 
the population size) and a cloud of asynchronous oscillators. The frequency of the cluster 
is larger than the frequency of the cloud oscillators. Noteworthy, these states appear in the 
case when autonomous oscillators are monostable, though close to the bistability domain (for 
sufficiently strong feedback, ar > 1, one time-delayed unit admits multiple periodic solu¬ 
tions). It means that the systems become bistable due to coupling and, thus, the chimeralike 
state emerges due to the dynamically sustained bistability. 
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VII. AN IMPORTANT APPLICATION: NEUROSCIENCE 


Within fifty years since Norbert Wiener hypothesized the role of the collective synchrony 
in the brain dynamics, importance of this phenomenon for neuroscience became highly recog¬ 
nized. It is now well-accepted that macroscopic neural rhythms appear due to coordination 
of bring and/or bursting of interconnected neuronal network and the Kuramoto model, in 
spite of its simplicity, became a paradigmatic and widely used setup in this held, see^ for 
a review. Noteworthy, the neuroscience applications turned out to be very fruitful for the 
theoretical development, posing quite interesting, from the viewpoint of nonlinear dynamics, 
problems. In particular, in the context of complex dynamical states which are of general 
interest, we have several times mentioned neuronal models, e.g. the Hindmarsh-Rose one. 
Below, without an ambition to provide a comprehensive review, we outline several relevant 
issues. 

For neuronal models, synchronization in a fully-connected network was observed for map- 
based and continuous time models, both for the regimes of periodic spiking and bursting, 
see e.g. Ref.—. A bit more detailed models consist of excitatory and inhibitory neuronal 
subpopulations^^. Another issue is a detailed description of synaptic coupling^, also with 
account of plasticity^. 

Fully connected network is certainly a rather crude approximation for a neuronal popu¬ 
lation. However, in some cases it is quite reasonable, as indicated by a numerical study of 
randomly coupled networks of map-based neurons^: if every unit is connected only to 0.5% 
of elements in the ensemble, then synchronization properties are practically indistinguishable 
from the fully connected case. For many problems the assumption of random connectivity is 
not appropriate; in this case one has to use physiologically motivated connectivity structure. 
An interesting approach has been recently suggested to treat large random networks of cou¬ 
pled oscillators^^. To adequately perform the thermodynamic limit and preserve disorder 
due to randomness of connections, a heterogeneous mean-beld approach has been developed 
in which disorder remains the same while the size of the system grows. This approach yields 
a description of both microscopic and global features of neuronal synchrony in the model^. 
Another interesting observation relates to diluted random networks of spike-coupled neu- 
rons2^: while for weak coupling synchrony establishes quite quickly, for large coupling a 
very long (in fact, exponential in the network size) transient disordered state is observed. 
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characterized by a negative largest Lyapunov exponent (so-called “stable chaos”). 


A. Control of collective synchrony 


As an example of a particular application of the discussed theoretical ideas we address 
the problem of suppression of the collective synchrony. This problem is relevant for a medi¬ 
cal technique, called deep brain stimulation (DBS). This technique is exploited to treat the 
Parkinson’s disease if it cannot be cured by medication. The DBS implies electrical stimula¬ 
tion of deep brain structures through the implanted micro-electrodes. The modern devices 
use the constant frequency (ca. 100 — 120Hz) pulse stimulation, which is typically applied 
around the clock. Although the exact mechanisms of DBS are not yet quite understood, it 
is widely used to reduce the limb tremor. 

Analysis of electrical or magnetic brain activity shows that Parkinsonian tremor is as¬ 
sociated with the pronounced spectral power at 10 to 12 Hence, it is reasonable to 

hypothesize that this pathological rhythm emerges due to synchronization in a neuronal pop¬ 
ulation and to consider the DBS as a desynchronization problent^. The main idea is then 
to develop efficient techniques which reliably suppress the unwanted activity with minimal 
stimulation. There are several directions in these studies. The hrst one implies open-loop 
stimulation with specially organized pulses^ through several closely spaced cites; the main 
assumption is that these pulses cause formation of several symmetrically positioned clusters, 
so that the mean held vanishes. 

Another direction in control of collective oscillation, based on the closed-loop feedback, 
was sug^sted in Refs. 94J]85| and verihed in experiments with electrochemical oscillators 
in Ref. j5|. The idea is quite transparent. Consider a globally coupled system: all elements 
are forced by the collective mode and synchronize due to this forcing. Suppose we measure 
the collective oscillation and feed it back, shifting its phase and properly amplifying it so 
that the feedback signal exactly compensates the mean held. In this case the oscillators 
become unforced and naturally desynchronize due to frequencies mismatch, internal noise, 
etc. Since the units desynchronize, the mean held tends to a constant, and so does the 
feedback signal. The constant component of the feedback signal can be easily eliminated; in 
this way one performs a vanishing stimulation control. It means that stimulation tends to 
zero as soon as the undesired oscillation is suppressed; this property is extremely important 
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for medical applications. The vanishing stimulation control can be implemented via a delayed 
feedback—!^ (see Ref. for a variant with nonlinear delayed feedback which is however 
not vanishing) or via a phase shifting passive system without delay^. Moreover, an adaptive 
strategy can be implementedSi, so that the desired state can be achieved in spite of unknown 
parameters of the system to be controlled. Notice that adjusting the phase shift introduced 
by the feedback loop one can, depending on the goal, suppress or enhance the collective 
synchrony. 


VIII. SURPRISES AND OUTLOOK 

In the hrst years after development of ideas of synchronization in large ensembles, es¬ 
pecially after Y. Kuramoto constructed his self-consistent theory, it looked like the very 
phenomenon of synchronization transition is rather simple and universal. In this paper we 
tried to show how in fact nontrivial is even the simplest setup: features like partial integra- 
bility, existence of an exact low-dimensional manifold, nontrivial transition scenaria even for 
unimodal distributions of frequencies, chimera states, self-organized quasiperiodicity occur 
already in the simplest case of sine-coupled phase oscillators. For generalizations of this 
model one observes a plethora of dynamical phenomena which is still far from being ex¬ 
hausted. Furthermore, addition of complexity to the basic model, e.g. by consideration of 
coupled oscillators on networks, results in further growth of the diversity of possible regimes. 

It seems to us that in the nearest future we will experience spreading of the synchro¬ 
nization theory far beyond its initial scope of nonlinear dissipative dynamical systems, e.g. 
to cover quantum objects^^. On the experimental level, advanced methods of oscillators’ 
control and of data analysis will possibly reveal microscopic details of nontrivial synchro¬ 
nization patterns. On the other hand, growing interest of mathematicians to these problems 
indicates that the held may become a part of mathematical physics as well. 
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